#include<cstdio>
#include<string.h>
#include<string>
#include<iostream>
#include<fstream>
#include<vector>
#include<cstdlib>
#include<set>
#include<map>
#include<stack>
#include<utility>
#include<climits>
#include<time.h>
#include<algorithm>
#include <boost/random.hpp>
#include <boost/random/uniform_01.hpp>
#include <boost/random/discrete_distribution.hpp>
#include"Graph.h"
using namespace std;
boost::mt19937 rng;
typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;


class Network{
public:
	Graph NodeToFactor;
	int nNode;
	vector<factor> Net;
	int *cardinality;

	Network(){
		Net = vector<factor>(); NodeToFactor = Graph(0); nNode = 0;
		cardinality = NULL;
	}
	Network(const vector<factor> net,const Graph G,const int N, const short *card){
		Net = net;
		NodeToFactor = G;
		nNode = N;
		cardinality = (int*)malloc(nNode*sizeof(int));
		for(int i=0;i<nNode;i++) cardinality[i]=card[i];
		//Get the initial degrees of each node

	}
//	Network& operator=(const Network &Netwk){
//		this->Net = Netwk.Net; this->NodeToFactor=Netwk.NodeToFactor; this->nNode = Netwk.nNode; this->cardinality = Netwk.cardinality;
//		return *this;
//	}
	~Network(){}
	//typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;

	void evidence(const vector<pair<int,short> > evid){

		for(unsigned i=0;i<evid.size();i++){
			vector<int> adjlist = NodeToFactor.getAdjacencyList(evid[i].first);
			float *condP;
			for (unsigned ii =0;ii<adjlist.size();ii++){
				factor F=Net[adjlist[ii]];
				factor T = factor();
				//Find the index of the element in the factor
				int index = find(F.first.begin(), F.first.end(), evid[i].first)-F.first.begin();
				//cout<<"Index: "<<index<<endl;
				if(index==0){
					//cout<<"Other Node:" <<F.first[1]<<endl;
					float tempfactor[cardinality[i]][cardinality[F.first[1]]];
						for(int ii=0;ii<cardinality[i];ii++){
							for(int jj= 0; jj<cardinality[F.first[1]];jj++ ){
								tempfactor[ii][jj]= F.second[ii*cardinality[F.first[1]] + jj].second;
							}
						}
					float sum = 0;
					condP = (float*)malloc(cardinality[F.first[1]]*sizeof(float));
					for(int k=0;k<cardinality[F.first[1]];k++)  {
						condP[k] = tempfactor[evid[i].second][k]; sum += condP[k];
					}
					for(int k=0;k<cardinality[F.first[1]];k++) condP[k]= condP[k]/sum;
					T.first.push_back(F.first[1]);
					for(int k=0;k<cardinality[F.first[1]];k++) {
						vector<short> temp = vector<short>(); temp.push_back(k);
						T.second.push_back(pair<vector<short>,float >(temp,condP[k]));
					}

				}else {
					//cout<<"Other Node:" <<F.first[0]<<endl;
					float tempfactor[cardinality[F.first[0]]][cardinality[i]];
						for(int ii= 0; ii<cardinality[F.first[0]];ii++ ){
							for(int jj=0;jj<cardinality[i];jj++){
								tempfactor[ii][jj]= F.second[ii*cardinality[i] + jj].second;
							}
						}
					float sum = 0;
					condP = (float*)malloc(cardinality[F.first[0]]*sizeof(float));
					for(int k=0;k<cardinality[F.first[0]];k++)  {
						condP[k]= tempfactor[k][evid[i].second]; sum += condP[k];
					}

					for(int k=0;k<cardinality[F.first[0]];k++) condP[k]= condP[k]/sum;
					T.first.push_back(F.first[0]);
					for(int k=0;k<cardinality[F.first[0]];k++) {
						vector<short> temp = vector<short>(); temp.push_back(k);
						T.second.push_back(pair<vector<short>,float >(temp,condP[k]));
					}


				}










				Net[adjlist[ii]]=T;


			}
			NodeToFactor.clear(evid[i].first);


		}

	}
	void Gibbs(const char *filename, const int nSample){
		//boost::random::discrete_distribution<> dist(P);
		cout<<"Inside Gibbs"<<endl;
		FILE *fsample = fopen(filename,"w+");
		cout<<fsample<<endl;
		//dist(rng);
		int *s = (int*)malloc(nNode*sizeof(int));

		int count = 0;

		for(int i=0;i<nNode;i++) { s[i]=2;  };
		//typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;

		for(int iSample =0;iSample<nSample; iSample++){


			for(int i=0;i<nNode;i++){
				vector<int> adjlist = NodeToFactor.getAdjacencyList(i);

				/************* Hard Code the Cardinality *****************/
				long double f[5]={0, 0, 0, 0, 0};
				for(unsigned j=0;j<adjlist.size();j++){
					int jngbr = adjlist[j];
					//showFactor(Net[jngbr]);

					if(Net[jngbr].first.size()==1){
						for(unsigned k=0;k<Net[jngbr].second.size();k++) f[k]*= Net[jngbr].second[k].second;
					}else{
						int index = find(Net[jngbr].first.begin(), Net[jngbr].first.end(), i)-Net[jngbr].first.begin();
						if(index==0){
							float tempfactor[cardinality[i]][cardinality[Net[jngbr].first[1]]];
							for(int ii=0;ii<cardinality[i];ii++){
								for(int jj= 0; jj<cardinality[Net[jngbr].first[1]];jj++ ){
									tempfactor[ii][jj]= Net[jngbr].second[ii*cardinality[Net[jngbr].first[1]] + jj].second;
								}
							}

//							for(int ii=0;ii<cardinality[i];ii++){
//								for(int jj= 0; jj<cardinality[Net[jngbr].first[1]];jj++ ){
//									cout<<tempfactor[ii][jj]<<" ";
//								}
//								cout<<endl;
//
//							}
							float *condP = (float*)malloc(cardinality[i]*sizeof(float));
							float sum = 0;
							for(int k=0;k<cardinality[i];k++)  {
								condP[k]= tempfactor[k][s[Net[jngbr].first[1]]] ; sum += condP[k];
							}
							assert(sum!=0);
							for(int k=0;k<cardinality[i];k++)  f[k] += log(condP[k]/sum);//f[k] *= condP[k]/sum;
							free(condP);




						}else {
							float tempfactor[cardinality[i]][cardinality[Net[jngbr].first[0]]];
							for(int ii= 0; ii<cardinality[Net[jngbr].first[0]];ii++ ){
								for(int jj=0;jj<cardinality[i];jj++){
									tempfactor[ii][jj]= Net[jngbr].second[ii*cardinality[i] + jj].second;
								}
							}

//							for(int ii=0;ii<cardinality[i];ii++){
//								for(int jj= 0; jj<cardinality[Net[jngbr].first[0]];jj++ ){
//									cout<<tempfactor[ii][jj]<<" ";
//								}
//								cout<<endl;
//
//							}
							float *condP = (float*)malloc(cardinality[i]*sizeof(float));
							float sum = 0;
							for(int k=0;k<cardinality[i];k++)  {
								condP[k]= tempfactor[s[Net[jngbr].first[0]]][k] ; sum += condP[k];
							}
							assert(sum!=0);
							for(int k=0;k<cardinality[i];k++)  f[k] += log(condP[k]/sum);//f[k] *= condP[k]/sum;
							free(condP);



						}


					}

				}

				long double sum = 0;
				for(int k=0;k<cardinality[i];k++)  {f[k]=exp(f[k]); sum +=f[k];}
				//cout <<sum <<endl;
				for(int k=0;k<cardinality[i];k++) f[k]=f[k]/sum;
				//cout<<"Node: "<<i<<" # Neighbours: "<<adjlist.size()<<endl;
				//for(int k=0;k<cardinality[i];k++) cout<<f[k]<<" ";



				boost::random::discrete_distribution<> dist(f);
				s[i] = dist(rng);
				//cout<<"Node: "<<i <<": ";
				if(i==0){for(int k=0;k<cardinality[i];k++) cout<<f[k]<<" ";
				cout <<endl;}
				//cout<<s[i]<<endl;

			}
			char *buff = (char*)malloc(sizeof(int));
			sprintf(buff,"%d",s[0]); string str =string(buff+'\0');
			for(int i=1;i<nNode;i++) {sprintf(buff,"%d",s[i]); str +=" " + string(buff+'\0'); }
			fprintf(fsample,"%s\n",str.c_str());




		}
		fclose(fsample);
	}
	factor getMarginal(const vector<int> elimorder,const int xNode){
		int count =0;
		for(unsigned i=0;i<elimorder.size();i++){
			if(xNode != elimorder[i] ){
				//if(count%1000==0) cout<<"Number of Nodes Deleted: "<<count<<endl;
				count++;
				if(!NodeToFactor.isDeleted[elimorder[i]]) BEL(elimorder[i]);
				//else{cout<<"ERROR: repeated node:"<<elimorder[i]<<endl; exit(EXIT_FAILURE);}
			}
		}
	    vector<int> factorlist = NodeToFactor.getAdjacencyList(xNode);

	    count =0;
		    //Get all the factors with the Node
		    //Multiply all the factors
		factor F = Net[factorlist[0]];
		for(unsigned i=1; i<factorlist.size();i++) F = product(F,Net[factorlist[i]]);

		float sum=0;
		for(unsigned j=0;j<F.second.size();j++) sum+=F.second[j].second;
		for(unsigned j=0;j<F.second.size();j++) F.second[j].second=F.second[j].second/sum;

		return F;

	}




	//typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;
	factor BEL(int xNode){
		//cout<<"******* Eliminating :"<<xNode<<"*******"<<endl;
	    vector<int> factorindices = NodeToFactor.getAdjacencyList(xNode);
	    //clock_t t= clock();

	    //Parallel Multiplication for Factors
	    vector<factor> factorlist = vector<factor>();

	    for(unsigned i=0;i<factorindices.size();i++) factorlist.push_back(Net[factorindices[i]]);
		vector<factor> temp1 = vector<factor>(); 	    vector<factor> temp2 = vector<factor>();

	    while(factorlist.size()!=1){
			//cout<<factorlist.size()<<endl;

			for(unsigned i=0; i<factorlist.size()/2;i++){
				temp1.push_back(factorlist[2*i]); temp2.push_back(factorlist[2*i+1]);
			}
			if(factorlist.size()%2==1) temp1.push_back(factorlist[factorlist.size()-1]);

			factorlist.clear();
			for(unsigned i=0;i<temp2.size();i++) factorlist.push_back(product(temp1[i],temp2[i]));
			if(temp1.size()>temp2.size()) factorlist.push_back(temp1[temp1.size()-1]);
			temp1.clear(); temp2.clear();
	    }
	    factor F = factorlist[0];

//	    factor F = Net[factorindices[0]];
//	    for(unsigned i=1; i<factorindices.size();i++) { F = product(F,Net[factorindices[i]]); }//showFactor(F);}
	    //clock_t t1= clock(); cout<<"Time for Product: " <<((float)(t1-t))/CLOCKS_PER_SEC <<endl;

	    //Sum over the node
	    factor S = sum(F,xNode);
	    NodeToFactor.clear(xNode);
	    int newFactorPosition = Net.size();
	    Net.push_back(S);
	    //clock_t t2= clock(); //cout<<"Time for Sum: " <<((float)(t2-t1))/CLOCKS_PER_SEC <<endl;


	    set<int> U;for(unsigned j=0;j<factorindices.size();j++) U.insert(factorindices[j]);

	    for(unsigned i=0; i<S.first.size();i++) {
		    //Link all remaining variables to the position of new factor
	    	NodeToFactor.insert(S.first[i],newFactorPosition);
	    	//cout<<S.first[i]<<": \n";
	    	//remove the edges from the remaining variables to the deleted factors
	    	vector<int>temp = NodeToFactor.getAdjacencyList(S.first[i]);
	    	for ( unsigned j=0;j<temp.size();j++) {
	    		if(U.find(temp[j])!=U.end()) {//cout<<S.first[i]<<"=>"<<it->second<<endl;
	    		NodeToFactor.remove(S.first[i],temp[j]); }
	    	}



	    }
	    U.clear();
	    //Put a null factor in place of all removed factor
	    vector<int> temp = vector<int>(); vector<pair<vector<short>,float > > tempvalueprob = vector<pair<vector<short>,float > >();
	    factor Empty(temp,tempvalueprob);
	    for(unsigned i=0;i<factorindices.size();i++) Net[factorindices[i]]=Empty;

	    //clock_t t3= clock(); //cout<<"Time for Elimination: " <<((float)(t3-t2))/CLOCKS_PER_SEC <<endl;

	    //showFactor(S);
	    return S;
	}

	factor product(factor f1, factor f2){
		//Find the Union of the two set of nodes
		set<int> S; set<int> U; //S is union, U is intersection of variables in two factors
		for(unsigned i=0;i<f1.first.size();i++) S.insert(f1.first[i]);
		for(unsigned i=0;i<f2.first.size();i++)
			if(S.find(f2.first[i])!=S.end()) U.insert(f2.first[i]);
		for(unsigned i=0;i<f2.first.size();i++) S.insert(f2.first[i]);


//		cout<<"Common Nodes: "; map<vectorkey,float> M;
//		for(set<int>::iterator it = U.begin();it!=U.end();it++) cout<<*it<<" ";
//		cout<<endl;

		int *index1 = (int*)malloc(U.size()*sizeof(int));
		int *index2 = (int*)malloc(U.size()*sizeof(int));


		int count = 0;
		for(set<int>::iterator it = U.begin();it!=U.end();it++){
			index1[count] = find(f1.first.begin(), f1.first.end(), *it) - f1.first.begin();
			index2[count] = find(f2.first.begin(), f2.first.end(), *it) - f2.first.begin();
			count ++;

		}



		vector<int> newfactorNodes=vector<int>();
		vector<pair<vector<short>,float > >  newfactorvalueprob = vector<pair<vector<short>,float > >();
		//Insert the nodes
		newfactorNodes = f1.first;
		for(unsigned i=0;i<f2.first.size();i++){
			if(U.find(f2.first[i])==U.end()){
				newfactorNodes.push_back(f2.first[i]);
			}
		}
		set<int> secondfactorIndex; count=0;
		for(set<int>::iterator it = U.begin();it!=U.end();it++){
			secondfactorIndex.insert(index2[count]); count ++;
		}
//typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;
		char *buff1 =(char*)malloc(sizeof(short)); char *buff2 =(char*)malloc(sizeof(short));

		for(unsigned i=0;i<f1.second.size();i++){
			for(unsigned j=0;j<f2.second.size();j++){
				int k = 0;string str1="", str2="";
				for(set<int>::iterator setit = U.begin();setit!=U.end();setit++){
					sprintf(buff1,"%u",f1.second[i].first[index1[k]]); str1+= string(buff1+'\0');
					sprintf(buff2,"%u",f2.second[j].first[index2[k]]); str2+= string(buff2+'\0');
					k++;
				}
				if(strcmp(str1.c_str(),str2.c_str())==0){
					//Copy the set of values from first pair
					vector<short> temp = f1.second[i].first;

					for(unsigned jj=0;jj<f2.second[j].first.size();jj++){
						if(secondfactorIndex.find(jj)==secondfactorIndex.end()){
							temp.push_back(f2.second[j].first[jj]);
						}
					}
					float prob = f1.second[i].second*f2.second[j].second;
					newfactorvalueprob.push_back(pair<vector<short>,float>(temp,prob));
				}



			}
		}

		factor F(newfactorNodes,newfactorvalueprob);
		free(index1); free(index2);
		return F;
	}
	//typedef pair<vector<int> ,vector<pair<vector<short>,float > > > factor;
	void BlockedGibbs(const vector<vector <int> > Trees){
		vector<Network> M = vector<Network>();
		for(unsigned i=0;i<Trees.size();i++){
			set <int> Nodes;
			for(unsigned j=0; j<Trees[i].size();j++){
				Nodes.insert(Trees[i][j]);
			}
			int nNodei = Trees[i].size();
			cout<<"Tree Size:"<<nNodei<<endl;
			short *cardi = (short*)malloc(nNodei*sizeof(short));
			for(unsigned j=0;j<Trees[i].size();j++) cardi[j] = cardinality[Trees[i][j]];
			Graph NTF= Graph(nNodei); vector<factor> Neti = vector<factor>();
			int count = 0;
			for(unsigned j=0;j<Trees[i].size();j++){
				//cout<<Trees[i][j]<<endl;
				vector<int> factorlist = NodeToFactor.getAdjacencyList(Trees[i][j]);
				for(unsigned k =0;k<factorlist.size();k++){
					if(this->Net[k].first.size()==1){
						Neti.push_back(Net[k]); NTF.insert(Trees[i][j],Neti.size()); count++;
					}else if(Nodes.find(this->Net[k].first[0])!=Nodes.end() && Nodes.find(this->Net[k].first[0])!=Nodes.end()){
						Neti.push_back(Net[k]); NTF.insert(Trees[i][j],Neti.size()); count++;
					}
				}
			}
			M.push_back(Network(Neti,NTF,nNodei,cardi));
			cout<<"Showing "<<i<<"th Tree: "<<endl;
			M[i].showDependency();
		}

	}
	factor sum(factor f,int xNode){
		map<vectorkey,float> M; vector<int> nodes=vector<int>();
		int index= find(f.first.begin(), f.first.end(), xNode) - f.first.begin();
		f.first.erase(f.first.begin()+index);

		for(unsigned i=0; i<f.second.size();i++){
			f.second[i].first.erase(f.second[i].first.begin()+index);
		}
		for(unsigned i=0; i<f.second.size();i++){
			if(M.find(f.second[i].first)==M.end()) M[f.second[i].first]=f.second[i].second;
			else M[f.second[i].first]+=f.second[i].second;
		}
		vector<pair<vector<short>,float > > value_prob = vector<pair<vector<short>,float > >();

		float sum=0;
		//Normalize the remaining values
		for(map<vectorkey,float>::iterator it=M.begin();it!=M.end();it++) sum+=it->second;
		for(map<vectorkey,float>::iterator it=M.begin();it!=M.end();it++){
			value_prob.push_back(pair<vector<short>,float >(it->first.key,it->second/sum));
		}



		return factor(f.first,value_prob);


	}

	void showFactor(factor F){
		cout<<"Connecting Node: ";
		for(unsigned j=0; j<F.first.size();j++){
			cout<<F.first[j]<<"  ";
		}
		cout<<endl<<"---------------"<<endl;
		for(unsigned j=0; j<F.second.size();j++){
			pair<vector<short>,float >  p = F.second[j];
			for(unsigned k=0;k<p.first.size();k++){
				cout<<p.first[k]<<"  ";
			}
		cout<<": "<<p.second<<endl;
		}
		cout<<endl<<endl;

	}
	void show(){
		cout<<"**********Printing the Factors**************"<<endl;

		for(unsigned i=0;i<Net.size();i++){
			cout<<"Factor: "<<i<<endl;
			cout<<"Connecting Node: ";
			for(unsigned j=0; j<Net[i].first.size();j++){
				cout<<Net[i].first[j]<<" ";
			}
			cout<<endl<<"---------------"<<endl;
			for(unsigned j=0; j<Net[i].second.size();j++){
				pair<vector<short>,float >  p = Net[i].second[j];
				for(unsigned k=0;k<p.first.size();k++){
					cout<<p.first[k]<<"\t";
				}
				cout<<": "<<p.second<<endl;
			}
			cout<<endl<<endl;
		}
	}
	void showDependency(){
		cout<<"**********Showing the Map in Between Nodes and Factors**************"<<endl;
		for(int i=0; i<nNode;i++){
		    vector<int> factorlist = NodeToFactor.getAdjacencyList(i);
		    cout << i << " =>";
		    for (unsigned j=0;j<factorlist.size();j++) cout << ' ' << factorlist[j];
		    std::cout << '\n';
		}

	}
};
